function [grdd,grdd3] =gradp_forward(f,x0,step,method)

    if nargin<4
        method='forward2';
    end

    f0 = feval(f,x0);
    [rf0,cf0] = size(f0);
    [rx0,cx0] = size(x0);
    grdd=zeros(rf0,rx0);

% Computation of stepsize (dh) for gradient 

    ax0 = abs(x0);
    if x0 ~= 0
        dax0 = x0./ax0;
    else
        dax0 = 1;
    end
%    dh = (1e-8)*max([ax0 (1e-2)*ones(rx0,1)]')';
    dh = (step).*max([ax0 (1e-2)*ones(rx0,1)]')';    
    dh = dh.*dax0;
    xdh = x0+dh;
    dh = xdh-x0;    % This increases precision slightly 
    arg = x0*ones(1,rx0);
    arg = arg-diag(diag(arg))+diag(xdh);
% Calculating f(xo+dh)          
    for i = 1:rx0
        disp(['Parameter # ' num2str(i) ' x0 ' num2str(x0(i)) ' x0+dh ' num2str(arg(i,i))])
        grdd(:,i) = feval(f,arg(:,i));
    
    end
    

    %% Forward 2 Method, default, quicker
    if lower(method)=='forward2'
        dh=repmat(dh',rf0,1);
        f0=repmat(f0,1,rx0);
        grdd = (grdd-f0)./dh;
    end
    
    %% Forward 3 Method, more precise
    if lower(method)=='forward3'
        xdh_ = x0+2*dh;
        arg_ = x0*ones(1,rx0);
        arg_ = arg-diag(diag(arg))+diag(xdh_);
% Calculating f(xo+2dh)   
        for i = 1:rx0
            disp(['Parameter # ' num2str(i) ' x0 ' num2str(x0(i)) ' x0+2dh ' num2str(arg_(i,i))])
                grdd_(:,i) = feval(f,arg_(:,i));
        end
    
    dh=repmat(dh',rf0,1);
    f0=repmat(f0,1,rx0);

    grdd = (4*grdd-grdd_-3*f0)./dh/2;
    end